#include <stdlib.h>
#include <stdio.h>

#include "rational.h"
#include "simplex.h"
#include "matrix.h"

/*
 *  Main Procedure Invoked for applying Simplex Method
 */
solution simplexProc(rational*** rationalMatrixWithArtVar, int rows, int cols, int** startingBasisVector, int* artificialVariableColumns, int numArtVar, rational** costVector) {
    
    rational** rationalMatrix;
    /* call function to create Rational Matrix 
     * and get artificial variables if any */
    solution phaseIResult, phaseIIResult;
    /* basis vector returned from phase I */
    /* int* phaseIBasisVector;
    getIntVector(&phaseIBasisVector, rows); */
    
    /* check if there are any artificial variables */
    
    if ( numArtVar > 0 ) {
    
        /* call phase I*/
        phaseIResult = phaseI(rationalMatrixWithArtVar, rows /*(rows+numArtVar)*/, cols, startingBasisVector, artificialVariableColumns, numArtVar );

        /* create New rational Matrix without art var */
        getRationalMatrix(&rationalMatrix, rows, (cols-numArtVar) ); /* does not contain the artificial variable columns */
        removeArtVarColumns(rationalMatrixWithArtVar, &rationalMatrix, artificialVariableColumns, numArtVar, rows, cols, startingBasisVector/*&phaseIBasisVector*/);
        
        /* replace given matrix with new matrix*/
        (*rationalMatrixWithArtVar) = rationalMatrix;
        
        printf("Matrix after Phase I :\n");
        printRationalMatrix(&rationalMatrix, rows, (cols-numArtVar));
        printf("Basis");
        printIntVector(startingBasisVector, rows);

        if ( phaseIResult == OPTIMAL ) {

            phaseIIResult = phaseII(&rationalMatrix, rows, (cols-numArtVar), startingBasisVector, costVector);

        } else {

            phaseIIResult = phaseIResult;
        }
    
    } else {
        
        /* proceed to phase II */
        phaseIIResult = phaseII(rationalMatrixWithArtVar, rows, cols, startingBasisVector, costVector);
        
    }
    
    return phaseIIResult;
    
}

/*
 * Phase-I Stage
 */
solution phaseI(rational*** rationalMatrix, int rows, int cols, int** startingBasisVector, int* artificialVariableColumns, int numArtVar) {
    /* artificialVariableColumns - column indices of artificial variables */
    int i;
    rational costValue;
    solution phaseIResult;
    /* create cost and basis vector for phase I */
    int *phaseIBasisVector;
    rational* phaseICostVector;
    getIntVector(&phaseIBasisVector, rows);
    getRationalVector(&phaseICostVector, cols); /* as long as rational Matrix */
    
    copyIntVector(startingBasisVector, &phaseIBasisVector, rows);
    
    for (i = 0; i < numArtVar; i++) {
        phaseICostVector[artificialVariableColumns[i]].numerator = -1;
        phaseICostVector[artificialVariableColumns[i]].denominator = 1;
    }
    
    printf("\nPHASE-I\n\n");
    
    /**/printf("Phase-I Cost");
    printRationalVector(&phaseICostVector, cols);
    
    phaseIResult = phaseIterations(rationalMatrix, rows, cols, &phaseIBasisVector, &phaseICostVector);
    
    if ( phaseIResult == OPTIMAL ) {
        /* check if artificial variables still present */
        /* costValue = basisCostVectorPdt(rows, &phaseIBasisVector, &phaseICostVector, rationalMatrix[cols]); */
        costValue = basisCostVectorPdt(rows, &phaseIBasisVector, &phaseICostVector, rationalMatrix, (cols-1));
        
        if (rE(costValue,itor(0))) {
            /* proceed to phase II */
            copyIntVector(&phaseIBasisVector, startingBasisVector, rows);
            return OPTIMAL;
        } else {
            /* infeasible */
            return INFEASIBLE;
        }
    }
    
}

/*
 * Phase-II Stage
 */
solution phaseII(rational*** rationalMatrix, int rows, int cols, int** startingBasisVector, rational** costVector) {
    
    int i;
    rational costValue;
    solution phaseIIResult;
    
    printf("\nPHASE-II\n\n");
    
    phaseIIResult = phaseIterations(rationalMatrix, rows, cols, startingBasisVector, costVector);
    
    return phaseIIResult;
    
}

/*
 * Procedure invoked to perform multiple iterations for getting more optimal solution
 * in Phase-I and Phase-II
 */
solution phaseIterations(rational*** rationalMatrix, int rows, int cols, int** startingBasisVector, rational** costVector) {
    
    int iteration = 0;
    solution iterationSolution = SUBOPTIMAL;
    /* create basis vectors for current and next iteration */
    int* currentBasisVector;
    getIntVector(&currentBasisVector, rows);
    copyIntVector(startingBasisVector, &currentBasisVector, rows);
    
    /* assumption : no. of iterations < MAX_ITERATIONS */
    while ( iterationSolution == SUBOPTIMAL ) {
        
        iterationSolution = singleIteration(rationalMatrix, rows, cols, &currentBasisVector, costVector, &currentBasisVector );
                
    }
    
    /* replace starting Basis Vector with optimal basis vector */
    if ( iterationSolution == OPTIMAL ) {
        copyIntVector(&currentBasisVector, startingBasisVector, rows);
    }    
    
    return iterationSolution;
}

/*
 * Procedure containing steps for single Iteration
 */
solution singleIteration(rational*** rationalMatrix, int rows, int cols, int** basisVector, rational** costVector, int** newBasisVector) {
    
    int column,row,i,isBasis;
    int enteringColumnCount,unbounded,rowIndex,optimalEnteringColumnVectorIndex;
    rational zj;
    rational cj__zj,theta,optimal_theta_cj__zj,theta_cj__zj;
    
    /* create matrix to store values for theta calculation */
    rational** thetaMatrix;
    getRationalMatrix(&thetaMatrix,rows,cols);
    
    
    /*
     * Find Entering Column
     */
    
    /* create vector to store columns that can enter */
    enteringColumn* enteringColumnVector;
    getEnteringColumnVector(&enteringColumnVector,cols);
    
    enteringColumnCount = -1; /* if negative at end, then no columns to enter */
    
    /* start checking among columns for entering column */
    for (column = 0; column < (cols-1); column++) { /* as last column is b_i */
                
        /* check if it is part of basis */
        isBasis = 0;
        for (i = 0; i < rows; i++) {
            if (column == (*basisVector)[i]) {
                isBasis = 1;
                break;
            }
        }
        if (isBasis) {
            continue;
        }
        
        /* optimality criterion check*/
        zj = basisCostVectorPdt(rows,basisVector,costVector,rationalMatrix,column);
        cj__zj = rSub((*costVector)[column], zj);
        
        if ( rGt(cj__zj, itor(0) ) ) {
            /* column can enter basis */
            enteringColumnCount++;
            enteringColumnVector[enteringColumnCount].columnIndex = column;
            enteringColumnVector[enteringColumnCount].cj_zj = cj__zj;
            
            /* check for unbounded*/
            unbounded = 1;
            /* for theta */
            rowIndex = -1;
            theta = itor(-1);
            
            for (row = 0; row < rows; row++) {
                if ( rGt((*rationalMatrix)[column][row],itor(0)) ) {
                    unbounded = 0; /* not unbounded for given column */
                    
                    /* compute Yik / Yrk */
                    thetaMatrix[column][row] = rDiv((*rationalMatrix)[(cols-1)][row]/*b_i*/, (*rationalMatrix)[column][row]);
                    
                    if (( rLt(thetaMatrix[column][row], theta) ) || ( rE(theta, itor(-1)) )) {
                        theta = thetaMatrix[column][row];
                        rowIndex = row;
                    }
                    
                }
            }
            
            if (unbounded == 1) {
                /* leave loop and return unbounded error condition */
                break;
            } else {
                /* theta exists */
                /*thetaVector[column] = theta;*/
                enteringColumnVector[enteringColumnCount].theta = theta;
                enteringColumnVector[enteringColumnCount].leavingRowIndex = rowIndex;
            }
            
        }

        
    }
    
    if (unbounded == 1) {
        /* error code to be returned */
        return UNBOUNDED;
    } else if (enteringColumnCount < 0 ) {
        /* code for optimality */
        return OPTIMAL;
    } else {
        /* find best entering column */
        
        optimalEnteringColumnVectorIndex = 0;
        optimal_theta_cj__zj = itor(0);
        
        for (i = 0; i < (enteringColumnCount+1); i++) {
            
            theta_cj__zj = rMul(enteringColumnVector[enteringColumnCount].theta,enteringColumnVector[enteringColumnCount].cj_zj);
            
            /* check if column optimal */
            if (rLt(optimal_theta_cj__zj, theta_cj__zj)) {
                /* optimal_theta_cj__zj = theta_cj__zj */
                optimal_theta_cj__zj.numerator = theta_cj__zj.numerator;
                optimal_theta_cj__zj.denominator = theta_cj__zj.denominator;
                
                optimalEnteringColumnVectorIndex = i;
            }
        }
        
        /* use optimal column index */        
    }
    
    /*
     * Transform input Matrix
     */
    
    for(i=0;i<15;i++) { printf("-"); }
    printf("\nPresent");
    printRationalMatrix(rationalMatrix, rows, cols);
    /* printf("Entering Column = %d \n",enteringColumnVector[optimalEnteringColumnVectorIndex].columnIndex);
    printf("Leaving Row = %d \n",enteringColumnVector[optimalEnteringColumnVectorIndex].leavingRowIndex); */
    printf("Entering Column -> %d \n",(enteringColumnVector[optimalEnteringColumnVectorIndex].columnIndex+1));
    printf("Leaving Row -> %d \n",(enteringColumnVector[optimalEnteringColumnVectorIndex].leavingRowIndex+1));
    
    changeRationalMatrix(rationalMatrix, rows, cols, enteringColumnVector[optimalEnteringColumnVectorIndex].columnIndex, enteringColumnVector[optimalEnteringColumnVectorIndex].leavingRowIndex );

    printf("\nChanged");
    printRationalMatrix(rationalMatrix, rows, cols);
    
    /* add new basis */
    (*newBasisVector)[enteringColumnVector[optimalEnteringColumnVectorIndex].leavingRowIndex] = enteringColumnVector[optimalEnteringColumnVectorIndex].columnIndex;
    
    printf("Basis");
    printIntVector(newBasisVector, rows);
    /* for(i=0;i<15;i++) { printf("-"); } */
    printf("\n");
    
    return SUBOPTIMAL;
    
}

void changeRationalMatrix(rational*** rationalMatrix, int rows, int cols, int enteringColumn, int leavingRow) {
    
    int column,row;
    rational** tmpRatMatrix;
    
    /* create tmp Matrix */
    getRationalMatrix(&tmpRatMatrix, rows, cols);
    
    for (column = 0; column < cols; column++) {  /* include b_i in processing */
        for (row = 0; row < rows; row++) {
            
            if ( row == leavingRow ) {
                
                /* Yrj' = Yrj/Yrk */
                /* (*rationalMatrix) */tmpRatMatrix[column][row] = rDiv( (*rationalMatrix)[column][row], (*rationalMatrix)[enteringColumn][leavingRow] );
                
            } else {
                
                /* Yij' = Yij - (Yik/Yrk)*Yrj */
                /* (*rationalMatrix) */tmpRatMatrix[column][row] =rSub(   (*rationalMatrix)[column][row], 
                                                        rMul(
                                                                rDiv(   (*rationalMatrix)[enteringColumn][row],
                                                                        (*rationalMatrix)[enteringColumn][leavingRow]
                                                                ),
                                                                (*rationalMatrix)[column][leavingRow]
                                                        ) 
                                                );
            }
        }
    }
    
    /* copy tmp matrix to given matrix */
    copyRationalMatrix(&tmpRatMatrix, rationalMatrix, 0, 0, 0, 0, rows, cols);
    
    return;        
    
}

void removeArtVarColumns(rational*** rationalMatrixWithArtVar, rational*** rationalMatrix, int* artificialVariableColumns, int numArtVar, int rows, int cols, int** basisVector) {
    
    int column,j,flag,basisIndex,newMatrixColumnCount;
    
    /* new basis vector */
    int* newBasisVector;
    getIntVector(&newBasisVector, rows);
    
    newMatrixColumnCount = 0;
    /* */
    for (column = 0; column < cols/*(rows+numArtVar)*/; column++ ) { /* as value column to be added too */
        
        /* check column for artificial variable column
         * if so, don't add
         */
        flag = 0;
        for (j = 0; j < numArtVar; j++) {
            if ( column == artificialVariableColumns[j] ) {
                flag = 1;
            }
        }
        if ( flag == 1 ) {
            continue;
        }
        
        /* check for basis column
         * if so, get their new column index
         */
        flag = 0;
        for (j = 0; j < rows; j++) {
            if ( column == (*basisVector)[j] ) {
                flag = 1;
                basisIndex = j;                        
            }
        }
        if ( flag == 1 ) {
            newBasisVector[basisIndex] = column;
        }
        
        /* add column */
        /*copyRationalVector(rationalMatrixWithArtVar[column], rationalMatrix[newMatrixColumnCount], rows);*/
        copyRationalMatrix(rationalMatrixWithArtVar, rationalMatrix, 0, column, 0, newMatrixColumnCount, rows, 1);
        newMatrixColumnCount++;
        
    }
    
    /* copy new basis to given basis vector */
    *basisVector = newBasisVector;
    
}

rational basisCostVectorPdt(int rank, int** basisVector, rational** costVector, rational*** rationalMatrix, int basisColumn) {
    
    int i;
    rational product, sum = itor(0);
    /* take product with basis columns */
    for (i = 0; i < rank; i++) {
        product = rMul((*costVector)[(*basisVector)[i]], (*rationalMatrix)[basisColumn][i]);
        sum = rAdd(sum,product);
    }
    
    /* return product; */
    return sum;
    
}

void getEnteringColumnVector(enteringColumn** rationalVector, int size) {
    
    enteringColumn* tmpVector;
    
    /* allocate space for matrix */
    tmpVector = (enteringColumn*) calloc(size, sizeof (enteringColumn)); /* initialize with 0 */
    
    *rationalVector = tmpVector;
    
}
